!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 4-4 order Runge-Kutta integrate subroutine
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Creator  : B. G.
! Date	  : 2015-12-05
! Revised  :
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Note :
!   I. all angle in rad
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine RK44(uBA1_RK,vBA1_RK,delPsi_RK,psi_RK,uBA2_RK,vBA2_RK)
implicit none
real*8::uBA1_RK,vBA1_RK,delPsi_RK,psi_RK,uBA2_RK,vBA2_RK
real*8::g,aAA2,u1,v1,psi1
real*8::k1,k2,k3,k4,l1,l2,l3,l4

g=1.4
! 1
u1=uBA1_RK
v1=vBA1_RK
psi1=psi_RK
aAA2=u1*u1+v1*v1                          ! M_asterisk**2
aAA2=(g+1)/2-(g-1)*aAA2/2                 ! (a/aA)**2
k1=v1*delPsi_RK
l1=(-u1+aAA2*(u1+v1/tan(psi1))/(v1*v1-aAA2))*delPsi_RK
! 2
u1=uBA1_RK+0.5*k1
v1=vBA1_RK+0.5*l1
psi1=psi_RK+0.5*delPsi_RK
aAA2=(g+1)/2-(g-1)*(u1*u1+v1*v1)/2
k2=v1*delPsi_RK
l2=(-u1+aAA2*(u1+v1/tan(psi1))/(v1*v1-aAA2))*delPsi_RK
! 3
u1=uBA1_RK+0.5*k2
v1=vBA1_RK+0.5*l2
aAA2=(g+1)/2-(g-1)*(u1*u1+v1*v1)/2
k3=v1*delPsi_RK
l3=(-u1+aAA2*(u1+v1/tan(psi1))/(v1*v1-aAA2))*delPsi_RK
! 4
u1=uBA1_RK+k3
v1=vBA1_RK+l3
psi1=psi_RK+delPsi_RK
aAA2=(g+1)/2-(g-1)*(u1*u1+v1*v1)/2
k4=v1*delPsi_RK
l4=(-u1+aAA2*(u1+v1/tan(psi1))/(v1*v1-aAA2))*delPsi_RK

uBA2_RK=uBA1_RK+(k1+2*k2+2*k3+k4)/6   
                                      ! radical non-d critical velocity along psi+delta_psi
vBA2_RK=vBA1_RK+(l1+2*l2+2*l3+l4)/6   
                                      ! circumferential non-d critical velocity along psi+delta_psi

end subroutine RK44
